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ABSTRACT 

The Milagro telescope monitors the northern sky for 100 GeV to 100 TeV transient emission through 
continuous very high energy wide-field observations. The large effective area and ^100 GeV energy 
threshold of Milagro allow it to detect very high energy (VHE) gamma-ray burst emission with much 
higher sensitivity than previous instruments, and a fluence sensitivity at VHE energies comparable to 
that of dedicated gamma-ray burst satellites at keV to MeV energies. Even in the absence of a positive 
detection, VHE observations can place important constraints on gamma-ray burst (GRB) progenitor 
and emission models. We present limits on the VHE flux of 40 s - 3 h duration transients nearby 
to earth, as well as sensitivity distributions which have been corrected for gamma-ray absorption by 
extragalactic background light and cosmological effects. The sensitivity distributions suggest that the 
typical intrinsic VHE fluence of GRBs is similar or weaker than the keV - MeV emission, and we 
demonstrate how these sensitivity distributions may be used to place observational constraints on the 
absolute VHE luminosity of gamma-ray bursts for any GRB emission and progenitor model. 
Subject headings: Gamma-Ray Bursts, Transient Search, TeV Gamma-Ray Observations 



1. INTRODUCTION 

Very high energy gamma-ray burst (GRB) observa- 
tions have the potential to constrain theoretical models 
of both the prompt and extended phases of GRB emis- 
sion. Models based on both internal and external shocks 
have predicted VHE fluence comparable to, or in cer- 
tain situations stronger than, the keV - MeV radiation, 
with durations ranging from shorter t han the keV - MeV 
burst to extended TeV afterglows iDer mer fc Chianel 
l2?XXl IPilla fc LoebllillMlZhang fc Mezarosll200lA . GRB 
emission above 100 GeV is particularly sensitive to the 
Lorentz factor and the photon density of the emitting 
material — and thus the distance of the radiating shock 
from the source — due to 77 — > e + e~ absorption in the 
emission region. 

EGRET observed several GRBs at GeV energies and 
observed no evi dence of a hig h energy rollover in the 
GRB spe c trum (jDingusI Il9 95). and recent results by 
!Gonzale3 l|2003D indicate that the spectrum of some 
GRBs contains a very hard, luminous, long duration 
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component. At energies above a few 100 GeV sev- 
eral observation attempts have been made, including: 
coincidence observa tions between the Tibet air-shower 
array and BATSE l)Amenomoril Il996|) . rapid follow-up 
observations by the Whi pple Air-Cherenkov telescope 
ijConnaughton et alJll99"^) . and coincident and monitor- 
ing studies by A IROBICC (jPadilla et alJll998|) . Whipple 
I Connaughton et all 119981). a nd the Milagro prototyp e 
1 Atkins et alJ l2000at iL^onorl 120001 iMcCulloughl [2001). 
While several of these studies have reported evidence for 
VHE GRB emission, no clear picture of the VHE emis- 
sion from gamma-ray bursts has emerged. 

In recent years the rate of possible GRB coincident 
detections has been relatively low due to the absence of 
sensitive, large field of view GRB detectors like BATSE. 
However, Milagro is a very sensitive detector with a low 
^100 GeV energy threshold and wide > 2 sr field of 
view, and is fully capable of autonomous identification 
of a VHE burst of emission. The low energy threshold 
is particularly important because of the reduced attenu- 
ation b y extragalactic backg round l ight near 100 GeV 
(iJellevI IT966t iPrimack et all Il999t iStecker fc de Jagerl 
1998), which dramatically increases the volume of space 
observed. This paper details the monitoring of the north- 
ern sky for VHE transients of 40 s - 3 hours duration with 
the Milagro observatory. 

2. THE 40 S TO 3 HOUR TRANSIENT SEARCH IN 
MILAGRO 

Very high energy gamma rays incident on the earth 
interact in the upper atmosphere and produce extensive 
air showers (EAS) which propagate to lower altitudes. 
Milagro uses the water Cherenkov technique to detect 
EASs by converting the wave front of relativistic parti- 
cles in the EAS into a front of Che renkov light, detecte d 
by photomultiplier tubes (PMTs) tl Atkins et al.ll2000b|) . 
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The amplitude and arrival time of the PMT signals are 
analyzed to determine the direction and characteristics 
of the ~ 1,800 EAS per second detected by Milagro, most 
of which are initiated by cosmic ray nuclei. The water 
Cherenkov technique enables Milagro to operate contin- 
uously. 

The precision of the Milagro event reconstruction de- 
pends strongly on the shower characteristics, leading to 
a point spread function (PSF) for the initiating particle 
direction which varies considerably from event to event. 
The arrival times measured by the PMTs are fit to a 
plane, after making corrections based on the pulse am- 
plitude in each PMT and the distance of the PMT from 
the reconstructed shower core. The shower fit proceeds 
through several iterations, at each stage removing times 
with large residuals to the fit shower plane. Showers 
reconstructed within 45° of zenith are kept for further 
analysis. The number of PMT times surviving in the fi- 
nal iteration is called n^t, which is used in conjunction 
with the reduced x 2 to characterize the quality of the 
reconstruction. The showers are then separated into 13 
classes, based on the values of and x 2 > sucn that 
events within each class have similar angular PSFs (see 
iMoralesI (J2002) for details). In general, the width of the 
PSF decreases as nst increases or x 2 decreases. The class 
with the worst PSF is removed from further analysis, at 
a negligible cost in sensitivity, because the PSF was so 
broad that it could not be readily handled by the code 
used to accumulate the sky maps. 

The characteristics of the light recorded by PMTs at 
the bottom of the reservoir can be used to reject back- 
ground cosmic ray showers in favor of gamma rays usin g 
a quantity we call "compactness" iAtkins et al.l [2003) . 
Events are assigned a weight, between and 1, based on 
the relative probability of the event being a background 
cosmic ray or a gamma ray, given the compactness value 
and to which of the 12 remaining PSF classes the event 
belongs. To speed the computation, events with weight 
less than 0.5 are rejected. 

This search imp l ement s an analysis method developed 
by IMorales et al.l l|2004|) which efficiently includes the 
event-by-event PSF and photon probability (i.e. back- 
ground rejection) information, enhancing the sensitivity 
of the observations. The technique uses the characteris- 
tics described above to determine a spatial photon prob- 
ability density — the expected PSF of the event multi- 
plied by the probability the shower was initiated by a 
gamma ray — which is then added to a sky map to form 
the estimated total photon density. Therefore, rather 
than considering each photon as a point at each recon- 
structed event location, the new technique instead places 
a smooth distribution at that location, with the shape of 
the distribution mirroring the reconstruction accuracy of 
the detector. The probability of the background, which 
primarily consists of cosmic ray events, producing an ob- 
served peak in the sky map is then given by the expected 
statistical distribution of the total photon density at that 
location on the sky map. The expected statistical distri- 
butions, which are typically not Gaussian in shape, can 
be calculated directly from the PSF characteristics of the 
detector and placed into a lookup table to enable very 
fast searches for unknown gamma-ray sources. This tech- 
nique is similar to maximum likelihood analyses in sen- 
sitivity, but is computationally fast enough to perform 



a real time transient search in a high data-rate detector 
like Milagro . 

The 40 s - 3 h transient search is performed by first 
dividing the data into 20 second intervals and creating 
sky maps, both for the actual data ("on source") and 
for the expected background in each time interval. From 
these initial maps, 40 s sky maps are formed by sum- 
ming all pairs of successive 20 s maps. The resulting 
set of sky maps is temporally oversampled by a factor 
of two, with every 20 s map used twice so that a 40 s 
burst is guaranteed to be largely contained in one of the 
overlapping 40 s maps. A search is performed on each of 
these maps for emission using the test statistic described 
below, then pairs of the 40 s maps are combined to make 
80 s maps (1/2 of the 40 s maps are discarded to main- 
tain the factor of two oversampling) , and so on. The 
nine specific time scales that were searched are 40, 80, 
160, 320, 640, 1280, 2560, 5120, and 10240 seconds. The 
total observation days for each time scale is shown in the 
second column of Table ^ The longer time scales have 
slightly less observation time due to gaps from occasional 
data dro pouts. Fu ll details of the transient search can be 
found in IMoralesI 1)20021) . 

The significance of excesses in the on-source sky maps, 
compared to the background sky maps, is evaluated using 
a test statistic. The test statistic is very closely related 
to the probability of an observation being produced by 
a background fluctuation, with approximations used to 
speed up the analysis breaking the exact one-to-one re- 
lationship and leading to a slight flattenin g of the test 
statistic distributions (Morales et al.ll2004|) . Since the 
value of the test statistic decreases with increasing tran- 
sient flux (lower probability of the observation being pro- 
duced by a background fluctuation), a discovery search 
simply looks for values of the test statistic below that 
expected from the number of independent observations. 
If a transient is identified, the statistical significance can 
be calculated from the test statistic. 

The data used for this analysis were taken between 
2001 May 2 nd and 2002 May 22 nd . Example distributions 
of the test statistic for 4 time scales are shown in Figure 
01 with no observed location for any of the nine search 
intervals having a test statistic below the ~10 -12 value 
expected from the number of independent observations. 
There is no evidence for transient VHE emission of 40 
s to 3 hours duration in the Milagro data between 2001 
May 2 nd and 2002 May 22 nd . 

3. LIMITS ON VHE GRB EMISSION 

Two sets of constraints are presented here — upper 
limits on the VHE flux at the earth, and sensitivity dis- 
tributions which constrain the absolute VHE luminosity 
of gamma-ray bursts. The flux limits have the advan- 
tage of being independent of cosmological models and 
directly comparable to previous VHE GRB limits. How- 
ever, due to the energy-dependent absorption of VHE 
photons by extragalactic background light (EBL), it is 
difficult to compare flux limits at the earth to theoreti- 
cal calculations. The sensitivity distributions place the 
observations into a theoretical context by correcting for 
absorption and cosmological effects. In conjunction with 
model predictions these can be used to determine the 
maximum VHE luminosity at the source which is consis- 
tent with the observations presented in this paper. 
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Fig. 1. — Test statistic distributions observed by Milagro between 
2001 May 2 nd and 2002 May 22 nd for 4 representative time scales. 
No evidence for transient VHE emission was observed. No test 
statistic below 10 — 12 was observed in any of the nine time scales 
studied. 



Since no observation in any of the nine time scales had 
a test statistic below 10 -12 , we can use this as a threshold 
to set a conservative limit on the maximum flux consis- 
tent with these observations. The relationship between 
the gamma ray flux and the test statistic was determined 
using a combination of Monte Carlo generated transient 
VHE signals and real data background. The simulated 
detector signal was modeled from 100 GeV to 21 TeV as 
a function of zenith angle and observed spectrum, prop- 
agated through the Milagro reconstruction code, added 
to representative data, and fed to the 40 s to 3 hour tran- 
sient detection program. If the test statistic for the sim- 
ulated signal fell below the 10 -12 threshold of this study, 
the signal would have been stronger than any observed 
fluctuation on any time scale. 

The flux limits represent the VHE transient signals 
which would have produced a signal with a test statistic 
below 10 -12 90% of the time (signal from a single burst 
which is excluded at the 90% confidence level). For this 
paper E~ 2,0 power-law spectra extending from 100 GeV 
to hard spectral cutoffs at 300 GeV, 1 TeV, or 21 TeV 
(as observed local to the earth) were chosen. These spec- 
tra serve as approximate models for an inverse Compton 
spectral bump at TeV energies or a hard VHE extension 
of the observed GRB s pectrum past the multi-GeV ob- 
servations of EGRET (Dingus 1995), with the spectral 
cutoff approximating either an intrinsic cutoff for nearby 
sources or EBL absorption for more distant sources. The 
limits are presented in Figure [2] as flux density as a func- 
tion of zenith angle and burst duration. 15 Given the 
~ 0.8 years of observation time and the ~ 2 sr field of 
view, these flux limits provide the strongest constraints 
to date on the the number of nearby bursts which emit 
TeV photons. For example, at a redshift of 0.024 (the 100 
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Fig. 2. — The 90% confidence upper limits for VHE transients of 
40 s duration and an E~ 2 ' spectrum for three different values of a 
hard spectral cutoff. The diamonds indicate the calculated limits 
on the normalization factor J in photons s cm" " 2 TeV -1 for the 



spectrum ^| = J( -^-y ) 2 U . Lines connecting the points have 
been drawn to guide the eye, with a solid line for a 21 TeV cutoff, 
dashed line for 1 TeV, and dash-dot line for 300 GeV. Sensitivity 
values for the eight longer time scales (80 s, 160 s, 320 s, 640 
s, 1280 s, 2560 s, 5120 s, and 10240 s) are not plotted, but the 
sensitivity improves as a J ex t~ °- 64 power law with the duration 
of the VHE transient. Monte Carlo statistics lead to an uncertainty 
in the upper limits of 19%. Systematic errors are due principally 
to uncertainties in the Monte Carlo simulation and are estimated 
to be +40%/-20%, for a total estimated error of +44%/-27%. 



15 The full set of limits for E~ 
cutoffs are shown in Morales! 1200: 
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Mpc distance of the nearby bursts identified by iNorris! 
(2002), with the cosmology used in that paper) our 90% 
confidence limit for a spectrum which extends to 21 TeV 
and burst emission which lasts 40 s at a 10.5 degree zenith 
angle is 1.1 x 10 46 ergs/s. 

However, at larger distances the usefulness of the limits 
from Figure |2 is complicated by the absorpt ion of VHE 
photons by extragalactic background light (lJellevlll966t 
iPrimack et aT]ll999l iStecker fc de Jagerlll99cf ~The ob- 
served spectrum is a convolution of the emitted spectrum 
and the extragalactic background light absorption — nei- 
ther of which are well determined. Limits on the absolute 
VHE luminosity of GRBs are thus necessarily model de- 
pendent. In setting the limits which follow, Milagro is 
more sensitive than might initially be expected because 
bursts above ~10 51 erg/s are sufficiently bright to be de- 
tected even after extinction of a few optical depths. 

In an effort to place the current observations in con- 
text, a set of assumptions about the emitted spectrum, 
EBL absorption, and cosmology have been chosen and 
sensitivity distributions calculated within this theoreti- 
cal framework. For a set of source luminosities and dura- 
tions, simulated GRBs were created with an E~ 2 emis- 
sion spectrum from 100 GeV to 21 TeV, an isotropic 
sky position, and following the star forma tion rate in z. 
The E BL absorpti on was determined by iBullock et alJ 
(|2004|) (similar to IPrimack et all |l999)) using recent 
results from semi-analytic modeling, and includes colli- 
sional starburst effec ts and a Kennicut initial mass func- 
tion l)Kennicutill983T) . 16 A ACDM cosmology (fi M = .3, 

16 The magn itude of the EBL absorptio n is quite uncertain, and 
the models by IStecker &: de Jagerl 119981) have considerably less 
attenuation. Consequently the limits presented here may be con- 
servative. 
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£1a = .7, h = 0.65) was used for both the EBL calculation 
and the cosmological corrections. The results are rather 
insensitive to the emitted spectrum and cutoff energy due 
to the strong EBL absorption which eliminates nearly all 
emission above 300 GeV at redshift « 0.3, and limits the 
distance to which Milagro can observe VHE emission to 
redshift < 0.7. The star formation rate (SFR) was mod- 
eled as lO 1 - 452 " 1 - 52 M©/Mpc 3 /year for z < 0.9 and 0.61 
M Q /Mpc 3 /year for z > 0.9 within each redsh ift bin, fol- 
lowing the discussion in lHopkins et alJ 1)200 lj) . 

The results of the sensitivity calculation are presented 
in TableQas the probability that Milagro would have ob- 
served a VHE GRB (i.e. test statistic < 10~ 12 ) as a func- 
tion of the absolute luminosity, 17 distance, and duration 
of the source. Due to increased atmospheric depth away 
from detector zenith, the low energy threshold of Milagro 
increases with zenith angle, as can be seen in Figure [21 
by the decreased sensitivity at large zenith angles when 
the emission is limited to lower energies. As GRBs be- 
come increasingly distant the EBL absorption effectively 
truncates the observed spectrum above a few hundred 
GeV. These effects lead to the effective field of view of 
Milagro becoming narrower with increased redshift, pro- 
ducing the gradual decline in detection probability seen 
in Tabled 

Given theoretical expectations for the GRB distance, 
luminosity, and duration distributions, the probabilities 
in Tabled can be used to calculate the upper limits pro- 
vided by these observations. For a 90% confidence up- 
per limit, the maximum GRB rate consistent with these 
observations is given by multiplying the expected num- 
ber of GRBs within each redshift bin by the detection 
probability in Table ^ an d adjusting the total number 
of GRBs until the mean number of events seen by Mila- 
gro is 2.3. If the model GRB distance distribution does 
not follow the SFR, the probabilities in Table Q] can still 
be used as approximate values, or corrections can be ap- 
plied to adjust for the distribution of events within each 
bin. The results are typically quoted as the number of 

17 The luminosity quoted is the isotropic luminosity assuming 
no beaming at the source, and is used to facilitate comparison with 
MeV to GeV measurements. 

18 One can use this fraction to roughly estimate that if all 
bursts were 10 51 erg/s, they would yield a flux brighter than 



GRBs/M Q of star formation or the average number of 
GRBs/Gpc 3 /year. 

As an example, consider a model which predicts that 
GRBs follow the star formation rate, with all GRBs 
emitting a characteristic 80 s pulse of VHE emission. 
The resulting 90% confidence upper limits for this model 
are 6.2 x 10~ 8 GRBs/M© of star formation (an average 
of 4.8 GRBs/ Gpc 3 / year over 0<z<0.5) for an isotropic 
luminosity of 10 51 ergs/s, or 1.1 x 10~ 8 GRBs/M Q of 
star formation (an average of 0.8 GRBs/Gpc 3 /year over 
0<z<0.5) for a luminosity of 10 52 ergs/s, with signifi- 
cantly tighter constraints if the GRB distribution trails 
the star formation rate (i.e. there are more low red- 
shift GRBs). Of the thirty-six GRBs with known dis- 
tances, five have a redshift below 0.5. 18 If the GRBs 
detected by BATSE follow the same distance distribu- 
tion, a rough estimate of the observed GRB rate yields 
-2.6 GRBs/Gpc 3 /year, or -3.4 x 10~ 8 GRBs/M Q of star 
formation if they follow the SFR. While detailed model 
calculations are needed to convert the probabilities in 
Table ^ into meaningful upper limits, comparison with 
the limits from this simple model suggests that if GRBs 
follow the star formation rate, the typical luminosity of 
40 s - 3 h VHE GRB counterparts is constrained to be 
similar to or less than the prompt keV - MeV emission. 
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the sensitivity in Figure 121 C~ 10~ 7 cm -2 s — - 1 ) out to redshift 
of about 0.5 for zenith angles less than ~ 30°, and approximately 
(5/36) X 0.84/(4vr) = 1% of all bursts would be detectible by Mi- 
lagro. 
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Table 1. Probability of GRB detection by Milagro. This table in- 
dicates THE PROBABILITY OF A GRB OF GIVEN DURATION, ABSOLUTE LUMI- 
nosity, and redshift within the 1.84 sr. field of view of this analysis 
being detected by the mllagro experiment. l x indicates the abso- 
lute luminosity (not corrected for beaming) in 10 x ergs/s between 
100 GeV and 21 TeV. Since the star formation rate and space volume 
both increase with redshift, the median redshift of the events used 
to calculate each probability is not in the center their redshift bin, 
but instead much nearer the upper redshift value. the triple dot 
symbols represent values which were not calculated due to computa- 
tional constraints. 
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